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Abstract. - We study energy propagation in locally time-periodically driven disordered nonlinear 
chains. For frequencies inside the band of linear Anderson modes, three different regimes are 
observed with increasing driver amplitude: 1) Below threshold, localized quasiperiodic oscillations 
and no spreading; 2) Three different regimes in time close to threshold, with almost regular 
oscillations initially, weak chaos and slow spreading for intermediate times, and finally strong 
diffusion; 3) Immediate spreading for strong driving. The thresholds are due to simple bifurcations, 
obtained analytically for a single oscillator, and numerically as turning-points of the nonlinear 
response manifold for a full chain. Generically, the threshold is nonzero also for infinite chains. 



Energy transport is rather well understood in lin- 
1 ear discrete Hamiltonan systems, where general solutions 
are linear combinations of eigenmodes, and the trans- 
port properties are completely determined by the na- 
ture of the spectrum of eigenfrequencies and the corre- 
sponding eigenmodes. For spatially periodic systems the 
linear spectrum is absolutely continuous and the eigen- 
modes are planewaves, implying that any initially localized 
wavepacket will generally disperse throughout the system, 
with amplitude vanishing to zero at infinite time. Con- 
versely, when the linear system is strongly disordered we 
have Anderson localization, characterized by a purely dis- 
crete linear spectrum with square summable (localized) 
eigenmodes. Then, initially localized wavepackets do not 
spread to zero and there is no energy diffusion. 

The aim of this paper is to understand how some as- 
pects of the transport properties of systems with linear 
Anderson localization are affected when nonlinearity is 
added. Obviously, for a nonlinear system, general solu- 
tions are not just linear combinations of eigenmodes. Still 
however, special (time-periodic) solutions generally exist, 
which continue single linear eigenmodes from the small- 
amplitude limit, and thus in some sense correspond to 



nonlinear eigenmodes. Due to the nonlinearities, the fre- 
quency of such a continued solution varies with ampli- 
tude, and each time it crosses a resonant linear eigenfre- 
quency a new peak appears at the location of the cor- 
responding linear mode. Since there are infinitely many 
resonant frequencies in any finite frequency interval, the 
continued solution develops infinitely many new peaks, 
and thus the strict continuation of a localized Anderson 
mode becomes a spatially extended mode for any non-zero 
amplitude [1]. Experimentally, such nonlinear derealiza- 
tion was reported in [2] for a two-dimensional fiber array. 
At larger amplitude, the spatial density of the peaks in- 
creases so that even in-between peaks the amplitude is 
never small, thus allowing the transportation of a sub- 
stantial energy current by phase torsion [3]. However, de- 
spite the existence of these extended solutions it has been 
shown [4-6] , that spatially localized time periodic solutions 
also persist with nonlinearities, but when the frequency is 
inside the linear band they are not in strict continuation 
of the linear modes. Instead, the frequency of a localized 
solution becomes a discontinuous function of its ampli- 
tude, with a gap located near each linear eigenfrequency. 
Although there are infinitely many gaps since the eigenfrc- 



p-1 



M. Johansson et al. 



qucncics form a dense set in the spectrum, the gap widths 
become exponentially small as functions of the spatial dis- 
tance of the resonance location, and as a consequence the 
"allowed" frequencies belong to a (fat) Cantor set with 
nonzero measure [4]. Indications of experimental excita- 
tion of a sequence of such "intraband discrete breathers" 
were given in [2] . Finally, when the nonlinearity is strong 
enough to drive the frequency outside the linear band, "ex- 
traband discrete breathers" (discrete solitons) are formed, 
and localization is generally enhanced as has been seen ex- 
perimentally for several types of photonic lattices [2, 7, 8] 

Thus, this indicates that nonlinearities on one hand 
maintain the existence of localized solutions which could 
trap energy forever, but on the other hand also gener- 
ate new extended solutions which could transport energy. 
To investigate the possible energy diffusion, numerical ex- 
periments have studied the time evolution of an initially 
localized wave packet [9-14]. This may be analyzed by 
expanding the solutions in the basis of Anderson modes, 
which then become coupled by nonlinear terms. It was 
early conjectured [9] that this coupling, if strong enough, 
would allow the energy to diffuse from one Anderson mode 
to another, leading to a (subdiffuse) spread throughout the 
system. Several numerical studies apparently confirmed 
that the second moment of the energy distribution of a sin- 
gle wave packet may diverge as a function of time [10, 11]. 
However, this does not necessarily imply that its ampli- 
tude vanishes at infinite time. By observing the time evo- 
lution of the participation number, which measures the 
localization length of the core of the wavepacket, it was 
shown [12] that for large enough initial amplitude the par- 
ticipation number remains finite, although the second mo- 
ment may diverge if some part of the energy spreads to 
infinity. The observations in [12] also suggested that the 
wavepacket may have a limit profile (possibly with infi- 
nite second moment) which is an almost periodic solution 
(with a discrete Fourier transform) and zero Lyapunov co- 
efficients analogous to KAM tori in finite systems. This 
conjecture is supported by rigorous results proving the ex- 
istence of such solutions in random DNLS systems [15,16]. 

This paper is devoted to another test for the energy 
transport property: the transmission through a system 
with a time periodic driving force locally applied. While 
the above arguments referred to finite energy wavepackcts 
in systems at zero temperature, the external driving force 
now acts as an energy source. We have investigated sev- 
eral nonlinear models with linear Anderson localization, 
but focus on two prototype examples: the Discrete Non- 
linear Schrodinger (DNLS) and Fermi-Pasta-Ulam (FPU) 
models, in both cases with driving at a single edge site 
of a semi-infinite chain. Numerically, we choose the chain 
length N to be sufficiently large to have no observable 
change when increasing N. Since we aim at describing a 
qualitative scenario generically valid for any typical disor- 
dered chain, we do not average over different realizations. 

A DNLS chain with random on-site potential and non- 
parametric periodic (harmonic) driving at a boundary site 



is described by the dynamical equations 

Itpn ~ V n 1p n + C(V>„+1 + Vn-l) + l^nl fyn = , (1) 

where C is the coupling constant, and the amplitude A n = 
AS n o of the driving force at frequency w is non-zero only 
at site n = 0. We take V n to be random variables with 
uniform distribution in the interval [— V, +V]. This gives 
a linear spectrum (for undriven chain) between [— 2C — 
V, +2C + V}. The explicit time-dependence of the force 
in eq. (1) is removed by transforming to a frame rotating 
with frequency lo, ip n (t) — (f> n (t)e lult , yielding 

i(j) n - {uj + V n )4> n + C{(j) n+ i +0„_i) + |</>„| 2 </>„ = A n . (2) 

For a given initial condition </> n (0), there are three inde- 
pendent parameters: A/C (driving strength), V/C (disor- 
der strength), oj/C (driving frequency). (We may always 
put C = 1 by rescaling time and wavefunction.) Equation 
(2) can be derived from the Hamiltonian 

n = {(w + v n M n \ 2 - c(^„ +1 < + r n+ iM 

n=0 

-^„I 4 + a„^ + a;^|, (3) 

since i(j) n = J^-, i</>* = — J^-. Thus, even in presence 
of driving H is a conserved quantity, characterizing the 
initial condition. If we choose 4> n {Q) = and apply the 
driving instantaneously at t = 0, then H = 0. However, 
if we slowly increase the force from zero to its final value 
A during a transient time, then TL is not constant during 
the transient, and generally H^O after the transient. We 
measure the "amount of excitation" in the system at a 
given time from the total norm (or power in optics) 

N(t)=Y J \4>n\\ (4) 

n 

which is a conserved quantity only in absence of driving. 

If all parameters |A/C|, |Vb/C|, \u>/C\ are large, we can 
neglect the coupling to the rest of the lattice and only 
consider the single anharmonic oscillator at n = 0: 

i<po-{u + V )<i>o + \<i>o\ 2 <f>o=A. (5) 

The nature of stationary solutions (constant \(f>o\) to eq. 
(5) is wellknown (see, e.g., ref. [17]). Due to conservation 
of H, it is integrable. Expressing <p — r e ie " (r > 
and 6o real) yields fo = — As'm9 ; 9o + w + Vo — tq = 
— (A/ro) cos 6q. Thus, a stationary solution (ro = 0) in 
presence of driving (A ^ 0) can only exist if 6*o =0 or 
7r, corresponding to oscillations inphased resp. antiphased 
to the driving force. They are obtained by solving the 
cubic equation r$ — ((J + Vo)ro =F A = 0, where upper 
(lower) sign corresponds to 9 = (ir). When A = 0, the 
solutions are ro = ("linear solution", H = 0) and ro = 
+ V if cj + Vq > (" nonlinear solution" , H = (lo + 
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Vb) 2 /2). For A > 0, the cubic equation for 9 = always 
has one (stable) solution for ro > 0, which continues the 
"nonlinear solution" towards larger r$ when u) + Vb > 
and the "linear solution" when u> + Vq < 0. In addition, 
for small A > and ui + Vo > the cubic equation for 
9 = 7r has two solutions with r > 0, one stable continuing 
the "linear solution" towards larger ro and one unstable 
continuing the "nonlinear solution" towards smaller ro. 
However, at a critical driving strength A = A^ = 2(u> + 
Vo) 3 / 2 / (3^3) these solutions bifurcate, so that for A > 
A(°' only the "nonlinear solution" with 9 = remains. 

The significance of this bifurcation is, that if A is in- 
creased slowly (adiabatically) from zero for a zero initial 
condition when uj + Vo > 0, the solution will follow the 
stable solution with 9 = n until it disappears at A = A^ c \ 
Above this threshold the solution cannot jump dynami- 
cally to the other stable stationary solution with 9 = 
since it has larger TL, but instead r becomes periodic with 
a new oscillation frequency (corresponding to a quasiperi- 
odic solution in the non-rotating frame). By contrast, for 
uj + Vb < there is always only one solution which contin- 
ues the "linear solution" for all A. 

For the full chain, an analogous threshold can be de- 
scribed by the nonlinear response manifold (NLRM) tech- 
nique, used in [5, 6] to calculate time-periodic breathers 
for a non-driven disordered lattice, and in [18] to ana- 
lyze transmission thresholds in driven nonrandom Klein- 
Gordon chains [19]. For the ID DNLS model, it can be 
simply implemented numerically. Looking for real (i.e., 
carrying no current) stationary solutions with <f> n = 0, eq. 
(2) can be written as a 2D map with 
which we may iterate backwards from n — N to n = 1, 



(uj + V n - x 2 n ) 
C 



Using as map initial (boundary) condition (xn,Vn) = 
(0, e) for a number of different e (typically ~ 1000 with 
max e ~ 10~ 12 for N ~ 50 — 100 sites) gives a number 
of points belonging to a manifold at the end of iteration 
at ((/>o , 4>\ ) . Combining this with the equation at the edge 
site n — 0, A — —(uj + Vb — Xq)xo + Cyo, and plotting, 
e.g., xo = 0o as a function of A gives a projection of the 
NLRM, from which we may obtain all real stationary so- 
lutions at given driving strength A as intersections of the 
NLRM with the vertical line at A. (The above set of map 
initial conditions yields only half the NLRM, the other 
half is obtained by adding (j) (—A) = — <j)o(A) but in order 
to not overload the figures it is not shown below) As dis- 
cussed in [18], each turning point (TP) of the NLRM can 
be associated with a threshold-like behaviour. Increasing 
the system size may yield more loops in the NLRM, but 
the structure of the first loops, and in particular the first 
TP, are generally found to remain unchanged for N > 100. 

We now illustrate numerically, how these NLRM TPs 
are related to transmission thresholds in the dynamics. 
Attempting to adiabatically follow the continuation of 
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Fig. 1: (a) (Half) NLRM for a random DNLS chain with 
V = 2.5 and u = 0. The apparent line through the origin 
is in fact a very thin loop, continuing the linear solution from 
the origin towards negative A until it reaches the first TP at 
(-1.9411, 0.697025) (inset in (a)). It then returns to the right, 
crosses A = at a small positive <f>o, reaches a second TP at 
(2.17, -0.97), and then escapes towards large negative values 
(not shown), (b) Time evolution of norm (main figure) and 
participation number (inset) for, from bottom to top, A = 1.5, 
A = 1.9 (just below first TP), A = 2.0 (between first and 
second TP), and A = 2.2 (above second TP). Transient time 
r = 100. (c), (d): Shape of the solution for A = 1.9, showing 
(c) the first transition from almost regular to slowly spreading, 
and (d) the second transition to strongly diffusive, behaviour. 



the linear stationary (i.e., with same frequency as driv- 
ing force) solution as long as it exists, we increase the 
driving strength slowly from zero to its final value, as 
A(t) = A(l - e-'/ T ), with typically r - 100 - 1000. We 
choose a rather strong disorder, V = 2.5, in order to have 
well localized linear modes (we put C = 1). To be con- 
crete, we show results obtained for one particular realiza- 
tion of the disorder, for which the on-site potential of the 
driving site is Vb ~ 2.48977. In fig. 1 we show (half) the 
NLRM and time-evolution of total norm (4) and participa- 
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Fig. 2: Adiabatic threshold A t h(o>), determined from the first 
TP of the NLRM, for the same DNLS chain as in Fig 1. The 
resolution in u> is 1CF 4 and chain length N = 100. 

tion number P = A/" 2 / |</>„| 4 for different A for a driving 
frequency in the middle of the linear band, u> = 0. 

When the driving strength is far below the first TP of 
the NLRM (A = 1.5 in fig. 1(b)), the response is essentially 
linear-like and, if the force is increased sufficiently slowly, 
the solution is essentially a stationary solution at the driv- 
ing frequency, localized close to the driving site and be- 
longing to the stable (lower) branch of the first (left) thin 
bifurcation loop in fig. 1(a)). The small (quasiperiodic) 
oscillations are due to non-adiabaticity, and will decrease 
if r is increased. Thus, this represents the sub-threshold- 
typc of non-transmitting dynamics. We have confirmed 
that this regular behaviour remains also for considerably 
longer integration times, e.g., 10 6 . 

For A just below the first TP (fig. 1(c), (d)), we observe 
a critical behaviour, where the small oscillations due to 
non-adiabaticity generally are weakly chaotic and cause 
the orbit to slowly diffuse in phase-space. At t <~ 5000 
it jumps to a different phase-space region, with bounded 
but larger oscillations (fig. 1(c)). We associate this with 
the existence of a neighboring stable stationary solution 
on the second loop of the NLRM (these solutions are close 
due to the thinness of the first loop). However, finally (t ~ 
13000, fig. 1(d)), the solution escapes into a continuously 
spreading, strongly chaotic state with rapid increase of Af 
and P. A similar scenario is observed also for A between 
the first and second TPs (A = 2.0 in fig. 1(b)), although 
in this case the weakly chaotic oscillations always develop 
as soon as A passes the first TP, even for larger r. 

Finally, for A above the second TP (A = 2.2 in fig. 1(b)), 
the spreading starts essentially immediately, and relevant 
measures like norm and participation number generally 
grow as t a , where a is close to 0.5 indicating a diffusive-like 
behaviour (for some different cases we obtained exponents 
between approximately 0.44 and 0.60). 

The main qualitative features discussed above are 
generic, although in particular the quantitative behaviour 
in the critical regime depends on the detailed structure 
of the NLRM, which for other parameter values may be 



considerably more complicated than in fig. 1(a). But gen- 
erally, for each random chain (fixed disorder realization 
and V) we can identify an adiabatic transmission thresh- 
old A t h(w) from the first TP of the NLRM, remembering 
that the real observed threshold is somewhat smaller due 
to nonadiabatic effects (finite r). The dependence A th (uj) 
exhibits sharp variations, as shown in fig. 2 for the chain 
used in fig. 1. For an infinite chain, Ath(w) is conjectured 
to be an upper semicontinuous function (larger than or 
equal to its limit at any point), which should become zero 
at any frequency resonant with a linear mode (i.e., on a 
dense set of points) , but nonzero for most frequencies and 
in average. In fig. 2 we can clearly observe the stronger 
resonances corresponding to modes located close to the 
driving site, while weaker resonances from distant modes 
should cause similar effects more difficult to observe nu- 
merically. We also confirmed that the transmission thresh- 
old decreases towards zero when approaching the lower 
band edge, as indicated by fig. 2. However, for frequen- 
cies smaller than the smallest resonance frequency (lower 
band edge for an infinite system, and u> w —2.7 for the 
finite system in fig. 2), there is no TP at all of the NLRM, 
and the numerical simulations show complete absence of 
transmission even for very large driving forces. 

To illustrate that the above scenario is indeed generic, 
and not due to special DNLS properties like conservation 
of H, we now discuss a random FPU chain with N sites 
and a cubic nonlinear force. The equations of motion are 

u n = K n (u n+ i - u n ) - K n -i(ll n - U n -\) 

+ - U n f - (U n -l - U n ) 3 , (6) 

where u n is the displacement of site n and K n are ran- 
dom coupling constants, independent and uniformly dis- 
tributed in [-?L m i n , if max]- In the following we fix if m ; n = 1 
and ifmax = 4 (we put Kq = 1). To simulate an im- 
pinging wave we impose on one end the boundary con- 
dition uo(t) = A cos Lut, while free boundary conditions, 
un+i = un, are enforced on the other side of the chain. 

An important difference between the FPU and DNLS 
models is that, with free boundary conditions on both 
sides, the FPU linear eigenmodes are well localized only 
for large enough frequencies, since the localization length 
£ diverges as £ ~ u~ 2 in the limit u> — > [20]. More pre- 
cisely, in a chain of N sites, all eigenvectors for which [20] 

/ 8(if-!) 
W < WL = )j N((K-*) - (K-i)*) (?) 

are in practice extended (phonon) modes. Therefore, we 
consider driving frequencies in the upper part of the linear 
spectrum, wl < w < w max where w ma x is the band-edge. 

In a series of numerical experiments we initialized the 
chain at rest (u n = u n = 0) and increased smoothly the 
driving amplitude from to A with a constant rate, uq = 
Acos(u)t) [l — e - '/ 7 "] , with typically r = 10. To simulate 
a semi-infinite lattice and observe a stationary state in 
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(U n — F(V n ,Un),U n ), where the function F is obtained 
by solving a cubic equation. In fig. 4 we show the last 
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Fig. 3: Average fluxes j versus amplitude driving A for different 
frequencies u. N = 512, Nd = 50, 7 = 0.1. The averages are 
performed on about 2 • 10 4 periods. Inset: same data in lin-log 
scale. Notice that the jump can be six orders of magnitude. 



the transmission regime, we steadily removed the energy 
injected from the driver by adding a damping term —"/u n 
to eq. (6) for a number Nd of rightmost sites (~ 10% of 
the total). As indicator, we used the local energy flux 



jn = \{ u n + U n+1 ) [K n (u n+1 



U n ) + («n+l _ u n) 3 



(8) 

After a transient of ~ 1000 driver periods, the site- and 
time-averaged energy flux j = (j2n=i* D in) / ( N ~ n d) 
was computed (the bar denotes time average). Figure 3 
shows results for several values of u (ui » wl ~ 0.44), 
for a fixed disorder realization. As in the DNLS case, 
a well-defined transmission threshold A t h exists for every 
frequency. Its exact value depends on the specific disorder 
realization (becoming small if there is a resonant mode 
close to the edge), but the qualitative behaviour is the 
same. Moreover, A t ^ is insensitive to variations of the size 
N (we checked sizes between N — 512 and N = 4096) and 
dissipation at damped sites Nd- Similarly to the DNLS 
case (cf. fig. 2 and the analytical result A^ for a single 
oscillator), there is an average tendency (disregarding the 
individual resonances) for A t h to increase with increasing 
frequency; note however that in the transmitting regime, 
the flux for a given driver amplitude increases with to. 

Also for FPU the threshold is related to a TP of the 
NLRM, which may be easily computed in a rotating- wave 
approximation (a more computationally expensive exact 
calculation could be done analogously to refs. [5,6,18]). 
Looking for solutions of the form u n (t) — U n cos uit and 
approximating cos 3 cut w (3/4) cos u>t in eq. (6), we get 

-uj 2 U n = K n (U n+1 - U n ) - K n _ x {U n - l/ n _i) 

+ \{U n +i-U n f -'^{U n ^-U n f. (9) 

This equation can be solved with respect to U n — U n -\ 
yielding a two dimensional backward map (J7 ra _i, V n -\) = 
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Fig. 4: Projection of the approximate FPU NLRM, obtained 
by iterating a few hundred points in an interval of size e ~ 10 -8 
around the origin, for to — 3.0237 and the same disorder real- 
ization as in fig. 3. Large points show the maximal oscillation 
amplitudes of the first particle, measured from simulation. 

iterate (Uo, Vq) (note that Uo = A) of a set of trajectories 
started along the unstable manifold of the origin. The 
curve is locally linear around A — (linear response) and 
then bends and turns wildly similar as discussed above. 
This computed manifold agrees very well with the data 
obtained from simulation when the driving is switched on 
very slowly (here r = 400). As can be seen, the first 
TP (TPi in fig. 4) is at A = 0.8756 which is in excellent 
agreement with the observed transmission threshold 0.878. 

Below threshold energy does not propagate and, as in 
the DNLS case, only a few sites close to the driven bound- 
ary oscillate. As shown in fig. 5 their motion is quasiperi- 
odic, with spectra displaying peaks at frequencies of the 
form mi fix + . . . + m p Q, p . Each site oscillates with a dif- 
ferent set of frequencies To characterize the localiza- 
tion properties of the quasiperiodic state we computed the 
time-averaged energy density e n , 



K 
2 



1 



[K+i(t 



Un) + V n ( 

)], (io) 



with V n (x) = K n x 2 /2 + x 4 /4. As seen in fig. 6, it displays 
a slow decay along the chain compatible with a power-law, 
e„ oc n~ 3 3 . A tail with similar exponent was found also 
when averaging over some different disorder realizations 
in the nontransmitting regime. Presumably, this is due to 
the presence of long- wavelength almost extended modes in 
the linear spectrum [20] ; in DNLS an exponential decay is 
generally seen for the sub-threshold quasiperiodic state. 

The interpretation of the threshold as a transition from 
quasiperiodicity to chaos can be seen in the Fourier spectra 
of u n as an immediate broadening of the lines for A > A t h- 
Asymptotically, for the transmitting state e„ is found to 
reach a given profile as in the ordered case [21], reminiscent 
of stationary heat transport with two thermal baths [22]. 
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Fig. 5: Fourier power spectra of the velocities u n at sites 
n = 1, 5, 10 for lo = 3.0237 and A = 0.850 (just below thresh- 
old). The insets show the Poincare sections (u n (kT),u n (kT)) 
(k integer) where T = 2n/uj is the driving period. 



In summary, we have shown the existence of well- 
defined transmission thresholds for generic classes of lo- 
cally time-periodically driven nonlinear disordered Hamil- 
tonian chains. An adiabatic threshold for the driver am- 
plitude can be defined from the associated NLRM, and is 
nonzero except at a discrete set of resonant frequencies. 
Applying the force nonadiabatically generates quasiperi- 
odic solutions in the nontransmitting regime, and lowers 
the transmission threshold. Beyond the threshold, energy 
is transmitted diffusively through a chaotic state. These 
thresholds could be directly observable, e.g., in optical 
waveguide arrays [8,23] and for Bose-Einstein condensates 
in disordered potentials [24] . Although we here focused on 
two particular models, we have found qualitatively simi- 
lar results also for random Klein-Gordon lattices [5,6, 18] 
and the parametrically driven DNLS model [25]. A re- 
cent work [26] also reported a transmission threshold in 
a model similar to eq. (1), but with dissipation at both 
edges. However, while ref. [26] argues that their thresh- 
old in average should decrease to zero in the limit of an 
infinite chain, a major conclusion from our work is that 
the threshold generically remains nonzero also in the limit 
N — > oo for any typical disorder realization. This follows 
also from the fact that the NLRM generally has a finite 
slope at the origin (except for the zero-measure set of res- 
onant frequencies), which was proven in ref. [6] under ba- 
sic smoothness assumptions. The existence of a smooth 
NLRM in the thermodynamic limit could possibly also be 
rigorously proven, at least for the DNLS case, as a conse- 
quence of the Rucllc-Oscledec theorem. 
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